pSERG CIRCADIAN DISTRIBUTION OF STATUS EPILEPTICUS

---

This analysis evaluates the circadian distribution of refractory status epilepticus episodes

Install packages needed for analysis

# install.packages("gmodels")
library(gmodels)

## Warning: package 'gmodels' was built under R version 3.4.4

# install.packages("gdata")
library(gdata)

## Warning: package 'gdata' was built under R version 3.4.4

## gdata: Unable to locate valid perl interpreter
## gdata: 
## gdata: read.xls() will be unable to read Excel XLS and XLSX files
## gdata: unless the 'perl=' argument is used to specify the location
## gdata: of a valid perl intrpreter.
## gdata: 
## gdata: (To avoid display of this message in the future, please
## gdata: ensure perl is installed and available on the executable
## gdata: search path.)

## gdata: Unable to load perl libaries needed by read.xls()
## gdata: to support 'XLX' (Excel 97-2004) files.

## 

## gdata: Unable to load perl libaries needed by read.xls()
## gdata: to support 'XLSX' (Excel 2007+) files.

## 

## gdata: Run the function 'installXLSXsupport()'
## gdata: to automatically download and install the perl
## gdata: libaries needed to support Excel XLS and XLSX formats.

## 
## Attaching package: 'gdata'

## The following object is masked from 'package:stats':
## 
##     nobs

## The following object is masked from 'package:utils':
## 
##     object.size

## The following object is masked from 'package:base':
## 
##     startsWith

# install.packages("plyr")
library(plyr)

## Warning: package 'plyr' was built under R version 3.4.4

# install.packages("devtools")
library(devtools)

## Warning: package 'devtools' was built under R version 3.4.3

# devtools::install_github('hadley/ggplot2')
library(ggplot2)

## Warning: package 'ggplot2' was built under R version 3.4.4

# install.packages("plotly")
library(plotly)

## Warning: package 'plotly' was built under R version 3.4.4

## 
## Attaching package: 'plotly'

## The following object is masked from 'package:ggplot2':
## 
##     last_plot

## The following objects are masked from 'package:plyr':
## 
##     arrange, mutate, rename, summarise

## The following object is masked from 'package:stats':
## 
##     filter

## The following object is masked from 'package:graphics':
## 
##     layout

# install.packages("cosinor")
library(cosinor)

## Warning: package 'cosinor' was built under R version 3.4.3

Load the database

## Load the data
pSERGcircadian <- read.csv("G:\\PROJECTS\\pSERG Circadian\\Data\\pSERG.csv")

Data cleaning and preparation of variables for analysis

## Keep patients only if they are a refractory SE case (not a control)
pSERGcircadian <- pSERGcircadian[which(pSERGcircadian$SE_GROUP == 'refractory_case'), ]
pSERGcircadian$SE_GROUP <- droplevels(pSERGcircadian$SE_GROUP)

## Transform date of status epilepticus into date format
pSERGcircadian$DATESEIZURE <- as.Date(pSERGcircadian$DATESEIZURE, format = "%m/%d/%Y")

## Order by date of status epilepticus
pSERGcircadian <- pSERGcircadian[order(pSERGcircadian$PATIENT_LABEL, pSERGcircadian$DATESEIZURE), ]

## Delete duplicate episodes from the same patient
pSERGcircadian <- pSERGcircadian[!duplicated(pSERGcircadian$PATIENT_LABEL), ]


#################Create variable hour of SE######################
## Transform time of hour to numeric and delete cases with unknown hour of SE
pSERGcircadian$TIMESEIZURE_HOURS <- as.numeric(as.character(pSERGcircadian$TIMESEIZURE_HOURS))

## Warning: NAs introduced by coercion

pSERGcircadian$TIMESEIZURE_HOURS <- ifelse(pSERGcircadian$TIMESEIZURE_HOURS>99, pSERGcircadian$TIMESEIZURE_HOURS%/%100, pSERGcircadian$TIMESEIZURE_HOURS)

## Delete patients with no numeric values for hour of seizure onset or no information
pSERGcircadian <- pSERGcircadian[!is.na(pSERGcircadian$TIMESEIZURE_HOURS), ]

## Rename variable
pSERGcircadian$hourofSE <- pSERGcircadian$TIMESEIZURE_HOURS

## Rename 0 to 24
pSERGcircadian$hourofSE[pSERGcircadian$hourofSE == 0] <- 24

## Exact hour
pSERGcircadian$TIMESEIZURE_MINUTES <- as.numeric(as.character(pSERGcircadian$TIMESEIZURE_MINUTES))

## Warning: NAs introduced by coercion

pSERGcircadian$time <- pSERGcircadian$TIMESEIZURE_HOURS + (pSERGcircadian$TIMESEIZURE_MINUTES / 60)


#######################Transform age into numeric and delete patients with unknown age####################
# Transform age into numeric
pSERGcircadian$G_T_STTS_PLPTCS_EPISODE_MONTHS <- as.numeric(as.character(pSERGcircadian$G_T_STTS_PLPTCS_EPISODE_MONTHS))

## Warning: NAs introduced by coercion

pSERGcircadian$G_T_STTS_PLPTCUS_EPISODE_YEARS <- as.numeric(as.character(pSERGcircadian$G_T_STTS_PLPTCUS_EPISODE_YEARS))
pSERGcircadian$ageyears <- pSERGcircadian$G_T_STTS_PLPTCUS_EPISODE_YEARS + (pSERGcircadian$G_T_STTS_PLPTCS_EPISODE_MONTHS/12)
pSERGcircadian <- pSERGcircadian[complete.cases(pSERGcircadian[ , "ageyears"]), ]


##########################Delete patients with unknown sex################################################
pSERGcircadian <- pSERGcircadian[which(pSERGcircadian$SEX == "male" | pSERGcircadian$SEX == "female"), ]


#######################Transform variables for analysis###############
## Create variable delay
pSERGcircadian$delay[grepl("delay", pSERGcircadian$PAST)] <- 1
pSERGcircadian$delay[!grepl("delay", pSERGcircadian$PAST)] <- 0

## Create variable palsy
pSERGcircadian$palsy[grepl("palsy", pSERGcircadian$PAST)] <- 1
pSERGcircadian$palsy[!grepl("palsy", pSERGcircadian$PAST)] <- 0
 
## Create variable febrile
pSERGcircadian$febrile[grepl("febrile", pSERGcircadian$PAST)] <- 1
pSERGcircadian$febrile[!grepl("febrile", pSERGcircadian$PAST)] <- 0

## Create variable prior epilepsy
pSERGcircadian$priorepilepsy[grepl("epi", pSERGcircadian$PAST)] <- 1
pSERGcircadian$priorepilepsy[!grepl("epi", pSERGcircadian$PAST)] <- 0

## Create variable status
pSERGcircadian$status[grepl("status", pSERGcircadian$PAST)] <- 1
pSERGcircadian$status[!grepl("status", pSERGcircadian$PAST)] <- 0

## Create variable none
pSERGcircadian$none[grepl("none", pSERGcircadian$PAST)] <- 1
pSERGcircadian$none[!grepl("none", pSERGcircadian$PAST)] <- 0

Demographics

############################DEMOGRAPHICS#######################
## Gender
CrossTable(pSERGcircadian$SEX)

## 
##  
##    Cell Contents
## |-------------------------|
## |                       N |
## |         N / Table Total |
## |-------------------------|
## 
##  
## Total Observations in Table:  300 
## 
##  
##           |    female |      male | 
##           |-----------|-----------|
##           |       134 |       166 | 
##           |     0.447 |     0.553 | 
##           |-----------|-----------|
## 
## 
## 
## 

## Age
nobs(pSERGcircadian$ageyears)

## [1] 300

summary(pSERGcircadian$ageyears)

##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
##  0.08333  1.28243  4.16667  5.99404  9.85417 20.74167

sd(pSERGcircadian$ageyears, na.rm = TRUE)

## [1] 5.337184

# Patients 18 year-old or older
pSERGcircadian$eighteenormore <- as.factor(ifelse(pSERGcircadian$ageyears >= 18, 1, 0))
CrossTable(pSERGcircadian$eighteenormore)

## 
##  
##    Cell Contents
## |-------------------------|
## |                       N |
## |         N / Table Total |
## |-------------------------|
## 
##  
## Total Observations in Table:  300 
## 
##  
##           |         0 |         1 | 
##           |-----------|-----------|
##           |       293 |         7 | 
##           |     0.977 |     0.023 | 
##           |-----------|-----------|
## 
## 
## 
## 

## Race
CrossTable(pSERGcircadian$RACE)

## 
##  
##    Cell Contents
## |-------------------------|
## |                       N |
## |         N / Table Total |
## |-------------------------|
## 
##  
## Total Observations in Table:  300 
## 
##  
##                                     |                              arabic |                               asian |           black_or_african_american | native_hawaiian_or_pacific_islander |                        not_reported | 
##                                     |-------------------------------------|-------------------------------------|-------------------------------------|-------------------------------------|-------------------------------------|
##                                     |                                  10 |                                  10 |                                  63 |                                   1 |                                  13 | 
##                                     |                               0.033 |                               0.033 |                               0.210 |                               0.003 |                               0.043 | 
##                                     |-------------------------------------|-------------------------------------|-------------------------------------|-------------------------------------|-------------------------------------|
## 
## 
##                                     |                             unknown |                               white | 
##                                     |-------------------------------------|-------------------------------------|
##                                     |                                  17 |                                 186 | 
##                                     |                               0.057 |                               0.620 | 
##                                     |-------------------------------------|-------------------------------------|
## 
## 
## 
## 

## Ethnicity
CrossTable(pSERGcircadian$ETHNICITY)

## 
##  
##    Cell Contents
## |-------------------------|
## |                       N |
## |         N / Table Total |
## |-------------------------|
## 
##  
## Total Observations in Table:  300 
## 
##  
##                        |     hispanic_or_latino | not_hispanic_or_latino |           not_reported |                unknown | 
##                        |------------------------|------------------------|------------------------|------------------------|
##                        |                     53 |                    218 |                     19 |                     10 | 
##                        |                  0.177 |                  0.727 |                  0.063 |                  0.033 | 
##                        |------------------------|------------------------|------------------------|------------------------|
## 
## 
## 
## 

## Delay
CrossTable(pSERGcircadian$delay)

## 
##  
##    Cell Contents
## |-------------------------|
## |                       N |
## |         N / Table Total |
## |-------------------------|
## 
##  
## Total Observations in Table:  300 
## 
##  
##           |         0 |         1 | 
##           |-----------|-----------|
##           |       146 |       154 | 
##           |     0.487 |     0.513 | 
##           |-----------|-----------|
## 
## 
## 
## 

## Palsy
CrossTable(pSERGcircadian$palsy)

## 
##  
##    Cell Contents
## |-------------------------|
## |                       N |
## |         N / Table Total |
## |-------------------------|
## 
##  
## Total Observations in Table:  300 
## 
##  
##           |         0 |         1 | 
##           |-----------|-----------|
##           |       270 |        30 | 
##           |     0.900 |     0.100 | 
##           |-----------|-----------|
## 
## 
## 
## 

## Febrile
CrossTable(pSERGcircadian$febrile)

## 
##  
##    Cell Contents
## |-------------------------|
## |                       N |
## |         N / Table Total |
## |-------------------------|
## 
##  
## Total Observations in Table:  300 
## 
##  
##           |         0 |         1 | 
##           |-----------|-----------|
##           |       269 |        31 | 
##           |     0.897 |     0.103 | 
##           |-----------|-----------|
## 
## 
## 
## 

## Prior epilepsy
CrossTable(pSERGcircadian$priorepilepsy)

## 
##  
##    Cell Contents
## |-------------------------|
## |                       N |
## |         N / Table Total |
## |-------------------------|
## 
##  
## Total Observations in Table:  300 
## 
##  
##           |         0 |         1 | 
##           |-----------|-----------|
##           |       152 |       148 | 
##           |     0.507 |     0.493 | 
##           |-----------|-----------|
## 
## 
## 
## 

## Status
CrossTable(pSERGcircadian$status)

## 
##  
##    Cell Contents
## |-------------------------|
## |                       N |
## |         N / Table Total |
## |-------------------------|
## 
##  
## Total Observations in Table:  300 
## 
##  
##           |         0 |         1 | 
##           |-----------|-----------|
##           |       239 |        61 | 
##           |     0.797 |     0.203 | 
##           |-----------|-----------|
## 
## 
## 
## 

## None
CrossTable(pSERGcircadian$none)

## 
##  
##    Cell Contents
## |-------------------------|
## |                       N |
## |         N / Table Total |
## |-------------------------|
## 
##  
## Total Observations in Table:  300 
## 
##  
##           |         0 |         1 | 
##           |-----------|-----------|
##           |       202 |        98 | 
##           |     0.673 |     0.327 | 
##           |-----------|-----------|
## 
## 
## 
## 

## Transform CONVULSIVEDURATION into numeric
pSERGcircadian$CONVULSIVEDURATION <- as.numeric(as.character(pSERGcircadian$CONVULSIVEDURATION))

## Warning: NAs introduced by coercion

pSERGcircadian$CONVULSIVEmin <- pSERGcircadian$CONVULSIVEDURATION
pSERGcircadian$CONVULSIVEhr <- pSERGcircadian$CONVULSIVEDURATION*60
pSERGcircadian$convulsivedurationmin <- ifelse(pSERGcircadian$CONVULSIVEDURATIONUNITS == "min", pSERGcircadian$CONVULSIVEmin, pSERGcircadian$CONVULSIVEhr)

## Duration of SE
nobs(pSERGcircadian$convulsivedurationmin)

## [1] 292

summary(pSERGcircadian$convulsivedurationmin)

##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max.     NA's 
##      0.0     60.0    120.0   2194.9    260.8 172800.0        8

sd(pSERGcircadian$convulsivedurationmin, na.rm = TRUE)

## [1] 13500.34

## Hospital onset
CrossTable(pSERGcircadian$HOSPITALONSET)

## 
##  
##    Cell Contents
## |-------------------------|
## |                       N |
## |         N / Table Total |
## |-------------------------|
## 
##  
## Total Observations in Table:  300 
## 
##  
##           |        no |       yes | 
##           |-----------|-----------|
##           |       202 |        98 | 
##           |     0.673 |     0.327 | 
##           |-----------|-----------|
## 
## 
## 
## 

## Type of status epilepticus
CrossTable(pSERGcircadian$TYPESTATUS)

## 
##  
##    Cell Contents
## |-------------------------|
## |                       N |
## |         N / Table Total |
## |-------------------------|
## 
##  
## Total Observations in Table:  300 
## 
##  
##              |   continuous | intermittent | 
##              |--------------|--------------|
##              |           94 |          206 | 
##              |        0.313 |        0.687 | 
##              |--------------|--------------|
## 
## 
## 
## 

Figures of the distribution of status epilepticus episodes

## Create dummy variable for events of SE
pSERGcircadian$events <- 1


################Every hour###############
## Draw the number of SE episodes in a plot with stacked units
distributionhour <- ggplot(pSERGcircadian, aes(x = hourofSE, y = events)) + geom_bar(stat = 'identity', color = "#0076c0", fill = "#67771a") + 
  scale_x_discrete(limits = c(1,24)) +
  theme(panel.background = element_rect(fill = "white"), 
        axis.text = element_text(size = 18, color = "black", face = "bold"),
        axis.title = element_text(size = 24, color = "black", face = "bold"),
        axis.ticks.length = unit(0.4, "cm"), axis.ticks = element_line(size = 1.2)) +
  scale_x_continuous(breaks = c(4, 8, 12, 16, 20, 24)) +
  scale_y_continuous(breaks = c(5, 10, 15, 20, 25)) +
  labs(x= "Time of day (hour)", y= "Number of rSE episodes")

## Scale for 'x' is already present. Adding another scale for 'x', which
## will replace the existing scale.

distributionhour

# Interactive graph
ggplotly(distributionhour)

################Every 2 hours###############
# Create variable every 2 hours
pSERGcircadian$groups2[pSERGcircadian$hourofSE == 1] <- "1-2"
pSERGcircadian$groups2[pSERGcircadian$hourofSE == 2] <- "1-2"

pSERGcircadian$groups2[pSERGcircadian$hourofSE == 3] <- "3-4"
pSERGcircadian$groups2[pSERGcircadian$hourofSE == 4] <- "3-4"

pSERGcircadian$groups2[pSERGcircadian$hourofSE == 5] <- "5-6"
pSERGcircadian$groups2[pSERGcircadian$hourofSE == 6] <- "5-6"

pSERGcircadian$groups2[pSERGcircadian$hourofSE == 7] <- "7-8"
pSERGcircadian$groups2[pSERGcircadian$hourofSE == 8] <- "7-8"

pSERGcircadian$groups2[pSERGcircadian$hourofSE == 9] <- "9-10"
pSERGcircadian$groups2[pSERGcircadian$hourofSE == 10] <- "9-10"

pSERGcircadian$groups2[pSERGcircadian$hourofSE == 11] <- "11-12"
pSERGcircadian$groups2[pSERGcircadian$hourofSE == 12] <- "11-12"

pSERGcircadian$groups2[pSERGcircadian$hourofSE == 13] <- "13-14"
pSERGcircadian$groups2[pSERGcircadian$hourofSE == 14] <- "13-14"

pSERGcircadian$groups2[pSERGcircadian$hourofSE == 15] <- "15-16"
pSERGcircadian$groups2[pSERGcircadian$hourofSE == 16] <- "15-16"

pSERGcircadian$groups2[pSERGcircadian$hourofSE == 17] <- "17-18"
pSERGcircadian$groups2[pSERGcircadian$hourofSE == 18] <- "17-18"

pSERGcircadian$groups2[pSERGcircadian$hourofSE == 19] <- "19-20"
pSERGcircadian$groups2[pSERGcircadian$hourofSE == 20] <- "19-20"

pSERGcircadian$groups2[pSERGcircadian$hourofSE == 21] <- "21-22"
pSERGcircadian$groups2[pSERGcircadian$hourofSE == 22] <- "21-22"

pSERGcircadian$groups2[pSERGcircadian$hourofSE == 23] <- "23-24"
pSERGcircadian$groups2[pSERGcircadian$hourofSE == 24] <- "23-24"

# Give this variable an order
pSERGcircadian$groups2 <- factor(pSERGcircadian$groups2, c("1-2", "3-4","5-6", "7-8","9-10", "11-12","13-14", "15-16","17-18", "19-20","21-22", "23-24"))

## Draw the number of SE episodes in a plot with stacked units
distribution2hours <- ggplot(pSERGcircadian[!is.na(pSERGcircadian$groups2), ], aes(x = groups2, y = events))+geom_bar(stat ='identity', color = "#0076c0", fill = "#67771a") + 
  theme(panel.background = element_rect(fill = "white"), 
        axis.text = element_text(size = 18, color = "black", face = "bold"),
        axis.title = element_text(size = 24, color = "black", face = "bold"),
        axis.ticks.length = unit(0.4, "cm"), axis.ticks = element_line(size = 1.2)) +
  scale_x_discrete() + scale_y_continuous(breaks = c(10, 20, 30, 40)) +
  labs(x= "Time of day (hour)", y= "Number of rSE episodes")
distribution2hours

# Interactive graph
ggplotly(distribution2hours)

################Every 3 hours###############
# Create variable every 3 hours
pSERGcircadian$groups3[pSERGcircadian$hourofSE == 1] <- "1-3"
pSERGcircadian$groups3[pSERGcircadian$hourofSE == 2] <- "1-3"
pSERGcircadian$groups3[pSERGcircadian$hourofSE == 3] <- "1-3"

pSERGcircadian$groups3[pSERGcircadian$hourofSE == 4] <- "4-6"
pSERGcircadian$groups3[pSERGcircadian$hourofSE == 5] <- "4-6"
pSERGcircadian$groups3[pSERGcircadian$hourofSE == 6] <- "4-6"

pSERGcircadian$groups3[pSERGcircadian$hourofSE == 7] <- "7-9"
pSERGcircadian$groups3[pSERGcircadian$hourofSE == 8] <- "7-9"
pSERGcircadian$groups3[pSERGcircadian$hourofSE == 9] <- "7-9"

pSERGcircadian$groups3[pSERGcircadian$hourofSE == 10] <- "10-12"
pSERGcircadian$groups3[pSERGcircadian$hourofSE == 11] <- "10-12"
pSERGcircadian$groups3[pSERGcircadian$hourofSE == 12] <- "10-12"

pSERGcircadian$groups3[pSERGcircadian$hourofSE == 13] <- "13-15"
pSERGcircadian$groups3[pSERGcircadian$hourofSE == 14] <- "13-15"
pSERGcircadian$groups3[pSERGcircadian$hourofSE == 15] <- "13-15"

pSERGcircadian$groups3[pSERGcircadian$hourofSE== 16] <- "16-18"
pSERGcircadian$groups3[pSERGcircadian$hourofSE== 17] <- "16-18"
pSERGcircadian$groups3[pSERGcircadian$hourofSE== 18] <- "16-18"

pSERGcircadian$groups3[pSERGcircadian$hourofSE== 19] <- "19-21"
pSERGcircadian$groups3[pSERGcircadian$hourofSE== 20] <- "19-21"
pSERGcircadian$groups3[pSERGcircadian$hourofSE== 21] <- "19-21"
 
pSERGcircadian$groups3[pSERGcircadian$hourofSE== 22] <- "22-24"
pSERGcircadian$groups3[pSERGcircadian$hourofSE== 23] <- "22-24"
pSERGcircadian$groups3[pSERGcircadian$hourofSE== 24] <- "22-24"

# Give this variable an order
pSERGcircadian$groups3 <- factor(pSERGcircadian$groups3, c("1-3", "4-6", "7-9", "10-12", "13-15", "16-18", "19-21", "22-24"))

## Draw the number of SE in a plot
distribution3hours <- ggplot(pSERGcircadian[!is.na(pSERGcircadian$groups3), ], aes(x = groups3, y = events)) + geom_bar(stat = 'identity', color = "#0076c0", fill = "#67771a") + 
  theme(panel.background=element_rect(fill = "white"), 
        axis.text = element_text(size = 18, color = "black", face = "bold"),
        axis.title = element_text(size = 24, color = "black", face = "bold"),
        axis.ticks.length = unit(0.4, "cm"), axis.ticks = element_line(size = 1.2)) +
  scale_x_discrete() + scale_y_continuous(breaks = c(10, 20, 30, 40, 50)) +
  labs(x = "Time of day (hour)", y = "Number of rSE episodes")
distribution3hours

# Interactive graph
ggplotly(distribution3hours)

################Every 4 hours###############
#Create variable every 4 hours
pSERGcircadian$groupsfour[pSERGcircadian$hourofSE == 1] <- "1-4"
pSERGcircadian$groupsfour[pSERGcircadian$hourofSE == 2] <- "1-4"
pSERGcircadian$groupsfour[pSERGcircadian$hourofSE == 3] <- "1-4"
pSERGcircadian$groupsfour[pSERGcircadian$hourofSE == 4] <- "1-4"

pSERGcircadian$groupsfour[pSERGcircadian$hourofSE == 5] <- "5-8"
pSERGcircadian$groupsfour[pSERGcircadian$hourofSE == 6] <- "5-8"
pSERGcircadian$groupsfour[pSERGcircadian$hourofSE == 7] <- "5-8"
pSERGcircadian$groupsfour[pSERGcircadian$hourofSE == 8] <- "5-8"

pSERGcircadian$groupsfour[pSERGcircadian$hourofSE == 9] <- "9-12"
pSERGcircadian$groupsfour[pSERGcircadian$hourofSE == 10] <- "9-12"
pSERGcircadian$groupsfour[pSERGcircadian$hourofSE == 11] <- "9-12"
pSERGcircadian$groupsfour[pSERGcircadian$hourofSE == 12] <- "9-12"

pSERGcircadian$groupsfour[pSERGcircadian$hourofSE == 13] <- "13-16"
pSERGcircadian$groupsfour[pSERGcircadian$hourofSE == 14] <- "13-16"
pSERGcircadian$groupsfour[pSERGcircadian$hourofSE == 15] <- "13-16"
pSERGcircadian$groupsfour[pSERGcircadian$hourofSE == 16] <- "13-16"

pSERGcircadian$groupsfour[pSERGcircadian$hourofSE == 17] <- "17-20"
pSERGcircadian$groupsfour[pSERGcircadian$hourofSE == 18] <- "17-20"
pSERGcircadian$groupsfour[pSERGcircadian$hourofSE == 19] <- "17-20"
pSERGcircadian$groupsfour[pSERGcircadian$hourofSE == 20] <- "17-20"

pSERGcircadian$groupsfour[pSERGcircadian$hourofSE == 21] <- "21-24"
pSERGcircadian$groupsfour[pSERGcircadian$hourofSE == 22] <- "21-24"
pSERGcircadian$groupsfour[pSERGcircadian$hourofSE == 23] <- "21-24"
pSERGcircadian$groupsfour[pSERGcircadian$hourofSE == 24] <- "21-24"

#Give this variable an order
pSERGcircadian$groupsfour <- factor(pSERGcircadian$groupsfour, c("1-4","5-8","9-12","13-16","17-20","21-24"))

##Draw the number of SE in a plot
distribution4hours <- ggplot(pSERGcircadian[!is.na(pSERGcircadian$groupsfour), ], aes(x = groupsfour, y = events)) + geom_bar(stat = 'identity', color = "#0076c0", fill = "#67771a") + 
  theme(panel.background = element_rect(fill = "white"), 
        axis.text = element_text(size = 18, color = "black", face = "bold"),
        axis.title = element_text(size = 24, color = "black", face = "bold"),
        axis.ticks.length = unit(0.4, "cm"), axis.ticks = element_line(size = 1.2)) +
  scale_x_discrete() + scale_y_continuous(breaks = c(10, 20, 30, 40, 50, 60, 70)) +
  labs(x = "Time of day (hour)", y = "Number of rSE episodes")
distribution4hours

# Interactive graph
ggplotly(distribution4hours)

Primary outcome

# Test (Kolmogorov-Smirnov) if the distribution of SE during the day is consistent with a uniform distribution
ks.test(pSERGcircadian$hourofSE, "punif", 0, 24)

## Warning in ks.test(pSERGcircadian$hourofSE, "punif", 0, 24): ties should
## not be present for the Kolmogorov-Smirnov test

## 
##  One-sample Kolmogorov-Smirnov test
## 
## data:  pSERGcircadian$hourofSE
## D = 0.085, p-value = 0.0262
## alternative hypothesis: two-sided

####################Distribution of SE episodes through the day###########################
# Distribution of SE episodes if they were evenly distributed
nobs(pSERGcircadian$events)

## [1] 300

nobs(pSERGcircadian$events)/24

## [1] 12.5

#####################################################################################
# Create another database with the circadian data
occurrences <- count(pSERGcircadian, "hourofSE")
# Make it a dataframe
circadiandata <- data.frame(occurrences)
circadiandata$numberSE <- circadiandata$freq
circadiandata <- circadiandata[ , c("hourofSE", "numberSE")]
#####################################################################################



# Cosinor analysis
SEfit <- cosinor.lm(numberSE ~ time(hourofSE), data = circadiandata, period = 24)
summary(SEfit)

## Raw model coefficients:
##             estimate standard.error lower.CI upper.CI p.value
## (Intercept)  12.5000         0.8293  10.8746  14.1254  0.0000
## rrr          -2.2350         1.1728  -4.5336   0.0636  0.0567
## sss           0.8876         1.1728  -1.4110   3.1862  0.4492
## 
## ***********************
## 
## Transformed coefficients:
##             estimate standard.error lower.CI upper.CI p.value
## (Intercept)  12.5000         0.8293  10.8746  14.1254  0.0000
## amp           2.4048         1.1728   0.1062   4.7034  0.0403
## acr          -0.3780         0.4877  -1.3339   0.5778  0.4382

# Draw the fitting curve 1h
## Draw the number of SE in a plot
circadianfigure <- ggplot.cosinor.lm(SEfit) + geom_bar(data = pSERGcircadian, aes(x = hourofSE, y = events), stat = 'identity', color = "#0076c0", fill = "#67771a", alpha = 0.8) + 
  scale_x_discrete(limits = c(1,24)) +
  theme(panel.background = element_rect(fill = "white"), 
        axis.text=element_text(size = 18, color = "black", face = "bold"),
        axis.title=element_text(size = 24, color = "black", face = "bold"),
        axis.ticks.length = unit(0.4, "cm"), axis.ticks = element_line(size = 1.2)) +
  scale_x_continuous(breaks = c(4, 8, 12, 16, 20, 24)) +
  scale_y_continuous(breaks = c(5, 10, 15, 20, 25)) +
  labs(x = "Time of day (hour)", y = "Number of rSE episodes") +
  geom_segment(aes(x = 0, xend = 24, y = SEfit$coefficients[1], yend = SEfit$coefficients[1])) +
  geom_segment(aes(x = 10.55, xend = 10.55, y = SEfit$coefficients[1], yend = (SEfit$coefficients[1] + SEfit$coefficients[2]))) +
  geom_segment(aes(x = 22.55, xend = 22.55, y = SEfit$coefficients[1], yend = (SEfit$coefficients[1] - SEfit$coefficients[2])))

## Scale for 'x' is already present. Adding another scale for 'x', which
## will replace the existing scale.

circadianfigure

# Interactive graph
ggplotly(circadianfigure)

#Draw the data and the fitting curve
circadianonly <- ggplot.cosinor.lm(SEfit) +
  geom_point(data = circadiandata, aes(hourofSE, numberSE), color = "#67771a", size = 4) +
  scale_x_discrete(limits = c(1,24)) + 
  theme(panel.background = element_rect(fill = "white"),
        axis.text = element_text(size = 18, color = "black", face = "bold"),
        axis.title = element_text(size = 24, color = "black", face = "bold"),
        axis.ticks.length = unit(0.4, "cm"), axis.ticks = element_line(size = 1.2),
        axis.line = element_line(color= "black", size = 1.2)) +
  scale_x_continuous(breaks = c(4, 8, 12, 16, 20, 24)) +
  labs(x = "Time of day (hour)", y = "Number of rSE episodes") +
  geom_segment(aes(x = 0, xend = 24, y = SEfit$coefficients[1], yend = SEfit$coefficients[1])) +
  geom_segment(aes(x = 10.55, xend = 10.55, y = SEfit$coefficients[1], yend = (SEfit$coefficients[1] + SEfit$coefficients[2]))) +
  geom_segment(aes(x = 22.55, xend = 22.55, y = SEfit$coefficients[1], yend = (SEfit$coefficients[1] - SEfit$coefficients[2])))

## Scale for 'x' is already present. Adding another scale for 'x', which
## will replace the existing scale.

circadianonly

# Interactive graph
ggplotly(circadianonly)

Comparison patients with and without a neurologic history

## Proportion of patients with some neurological history per time bin
neurologichistory1h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 1, ]$none)[1] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 1, ]$none)[1] + table(pSERGcircadian[pSERGcircadian$hourofSE == 1, ]$none)[2])
neurologichistory2h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 2, ]$none)[1] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 2, ]$none)[1] + table(pSERGcircadian[pSERGcircadian$hourofSE == 2, ]$none)[2])
neurologichistory3h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 3, ]$none)[1] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 3, ]$none)[1] + table(pSERGcircadian[pSERGcircadian$hourofSE == 3, ]$none)[2])
neurologichistory4h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 4, ]$none)[1] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 4, ]$none)[1] + table(pSERGcircadian[pSERGcircadian$hourofSE == 4, ]$none)[2])
neurologichistory5h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 5, ]$none)[1] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 5, ]$none)[1] + table(pSERGcircadian[pSERGcircadian$hourofSE == 5, ]$none)[2])
neurologichistory6h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 6, ]$none)[1] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 6, ]$none)[1] + table(pSERGcircadian[pSERGcircadian$hourofSE == 6, ]$none)[2])
neurologichistory7h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 7, ]$none)[1] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 7, ]$none)[1] + table(pSERGcircadian[pSERGcircadian$hourofSE == 7, ]$none)[2])
neurologichistory8h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 8, ]$none)[1] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 8, ]$none)[1] + table(pSERGcircadian[pSERGcircadian$hourofSE == 8, ]$none)[2])
neurologichistory9h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 9, ]$none)[1] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 9, ]$none)[1] + table(pSERGcircadian[pSERGcircadian$hourofSE == 9, ]$none)[2])
neurologichistory10h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 10, ]$none)[1] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 10, ]$none)[1] + table(pSERGcircadian[pSERGcircadian$hourofSE == 10, ]$none)[2])
neurologichistory11h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 11, ]$none)[1] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 11, ]$none)[1] + table(pSERGcircadian[pSERGcircadian$hourofSE == 11, ]$none)[2])
neurologichistory12h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 12, ]$none)[1] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 12, ]$none)[1] + table(pSERGcircadian[pSERGcircadian$hourofSE == 12, ]$none)[2])
neurologichistory13h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 13, ]$none)[1] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 13, ]$none)[1] + table(pSERGcircadian[pSERGcircadian$hourofSE == 13, ]$none)[2])
neurologichistory14h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 14, ]$none)[1] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 14, ]$none)[1] + table(pSERGcircadian[pSERGcircadian$hourofSE == 14, ]$none)[2])
neurologichistory15h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 15, ]$none)[1] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 15, ]$none)[1] + table(pSERGcircadian[pSERGcircadian$hourofSE == 15, ]$none)[2])
neurologichistory16h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 16, ]$none)[1] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 16, ]$none)[1] + table(pSERGcircadian[pSERGcircadian$hourofSE == 16, ]$none)[2])
neurologichistory17h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 17, ]$none)[1] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 17, ]$none)[1] + table(pSERGcircadian[pSERGcircadian$hourofSE == 17, ]$none)[2])
neurologichistory18h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 18, ]$none)[1] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 18, ]$none)[1] + table(pSERGcircadian[pSERGcircadian$hourofSE == 18, ]$none)[2])
neurologichistory19h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 19, ]$none)[1] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 19, ]$none)[1] + table(pSERGcircadian[pSERGcircadian$hourofSE == 19, ]$none)[2])
neurologichistory20h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 20, ]$none)[1] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 20, ]$none)[1] + table(pSERGcircadian[pSERGcircadian$hourofSE == 20, ]$none)[2])
neurologichistory21h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 21, ]$none)[1] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 21, ]$none)[1] + table(pSERGcircadian[pSERGcircadian$hourofSE == 21, ]$none)[2])
neurologichistory22h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 22, ]$none)[1] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 22, ]$none)[1] + table(pSERGcircadian[pSERGcircadian$hourofSE == 22, ]$none)[2])
neurologichistory23h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 23, ]$none)[1] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 23, ]$none)[1] + table(pSERGcircadian[pSERGcircadian$hourofSE == 23, ]$none)[2])
neurologichistory24h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 24, ]$none)[1] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 24, ]$none)[1] + table(pSERGcircadian[pSERGcircadian$hourofSE == 24, ]$none)[2])

circadiandata$neurologic <- c(neurologichistory1h, neurologichistory2h, neurologichistory3h, neurologichistory4h, neurologichistory5h, neurologichistory6h, neurologichistory7h, neurologichistory8h, neurologichistory9h, neurologichistory10h, neurologichistory11h, neurologichistory12h, neurologichistory13h, neurologichistory14h, neurologichistory15h, neurologichistory16h, neurologichistory17h, neurologichistory18h, neurologichistory19h, neurologichistory20h, neurologichistory21h, neurologichistory22h, neurologichistory23h, neurologichistory24h)
circadiandata$neurologic <- ifelse(is.na(circadiandata$neurologic) == TRUE, 0, circadiandata$neurologic)
circadiandata$numberneurologic <- circadiandata$numberSE * circadiandata$neurologic


## Cosinor analysis
SEfitcomparisonneurologic <- cosinor.lm(numberSE ~ time(hourofSE) + numberneurologic + amp.acro(numberneurologic), data = circadiandata, period = 24)
summary(SEfitcomparisonneurologic)

## Raw model coefficients:
##                      estimate standard.error lower.CI upper.CI p.value
## (Intercept)            3.3774         1.2213   0.9837   5.7710  0.0057
## numberneurologic       1.0820         0.1477   0.7926   1.3714  0.0000
## rrr                    0.1795         1.3358  -2.4387   2.7977  0.8931
## sss                   -0.2694         2.1349  -4.4537   3.9149  0.8996
## numberneurologic:rrr  -0.2406         0.1527  -0.5398   0.0587  0.1151
## numberneurologic:sss  -0.0598         0.2469  -0.5436   0.4240  0.8086
## 
## ***********************
## 
## Transformed coefficients:
##                            estimate standard.error lower.CI upper.CI
## (Intercept)                  3.3774         1.2213   0.9837   5.7710
## [numberneurologic = 1]       1.0820         0.1477   0.7926   1.3714
## amp                          0.3237         1.9249  -3.4489   4.0964
## [numberneurologic = 1]:amp   0.3348         1.8762  -3.3425   4.0121
## acr                         -0.9831         5.0164 -10.8151   8.8488
## [numberneurologic = 1]:acr   1.3874         3.6858  -5.8366   8.6113
##                            p.value
## (Intercept)                 0.0057
## [numberneurologic = 1]      0.0000
## amp                         0.8664
## [numberneurologic = 1]:amp  0.8584
## acr                         0.8446
## [numberneurologic = 1]:acr  0.7066

test_cosinor(SEfitcomparisonneurologic, "numberneurologic", param = "amp")

## Global test: 
## Statistic: 
## [1] 0
## 
## 
##  P-value: 
## [1] 0.9907
## 
##  Individual tests: 
## Statistic: 
## [1] 0.01
## 
## 
##  P-value: 
## [1] 0.9907
## 
##  Estimate and confidence interval[1] "0.01 (-1.86 to 1.88)"

test_cosinor(SEfitcomparisonneurologic, "numberneurologic", param = "acr")

## Global test: 
## Statistic: 
## [1] 0.25
## 
## 
##  P-value: 
## [1] 0.6163
## 
##  Individual tests: 
## Statistic: 
## [1] 0.5
## 
## 
##  P-value: 
## [1] 0.6163
## 
##  Estimate and confidence interval[1] "2.37 (-6.9 to 11.64)"

##Graph
neurologichistoryfigure <- ggplot.cosinor.lm(SEfitcomparisonneurologic, x_str = "numberneurologic") +
  geom_bar(data = pSERGcircadian[pSERGcircadian$none == 0, ], aes(x = hourofSE, y = events), stat = 'identity', color = "#5698a3", fill = "#5698a3", alpha = 0.2) + 
  geom_bar(data = pSERGcircadian[pSERGcircadian$none == 1, ], aes(x = hourofSE, y = events), stat = 'identity', color = "#a30234", fill = "#a30234", alpha = 0.2) +
    scale_x_discrete(limits = c(1,24)) +
  theme(panel.background=element_rect(fill = "white"), 
        axis.text=element_text(size = 18, color = "black", face = "bold"),
        axis.title=element_text(size = 24, color = "black", face = "bold"),
        axis.ticks.length=unit(0.4, "cm"), axis.ticks = element_line(size = 1.2),
        legend.position = "none") +
  scale_x_continuous(breaks = c(4, 8, 12, 16, 20, 24)) +
  scale_y_continuous(breaks = c(5, 10, 15, 20)) +
  labs(x = "Time of day (hour)", y = "Number of rSE episodes") +
  geom_segment(aes(x = 0, xend = 24, y = SEfitcomparisonneurologic$coefficients[1], yend = SEfitcomparisonneurologic$coefficients[1]), color = "#a30234") +
  geom_segment(aes(x = 20.2, xend = 20.2, y = SEfitcomparisonneurologic$coefficients[1], yend = (SEfitcomparisonneurologic$coefficients[1] + SEfitcomparisonneurologic$coefficients[3])), color = "#a30234") +
  geom_segment(aes(x = 8.2, xend = 8.2, y = SEfitcomparisonneurologic$coefficients[1], yend = (SEfitcomparisonneurologic$coefficients[1] - SEfitcomparisonneurologic$coefficients[3])), color = "#a30234") +
  geom_segment(aes(x = 0, xend = 24, y = SEfitcomparisonneurologic$coefficients[1] + SEfitcomparisonneurologic$coefficients[2], yend = SEfitcomparisonneurologic$coefficients[1] + SEfitcomparisonneurologic$coefficients[2]), color = "#5698a3") +
  geom_segment(aes(x = 17.25, xend = 17.25, y = SEfitcomparisonneurologic$coefficients[1] + SEfitcomparisonneurologic$coefficients[2], yend = (SEfitcomparisonneurologic$coefficients[1] + SEfitcomparisonneurologic$coefficients[2] + SEfitcomparisonneurologic$coefficients[4])), color = "#5698a3") + geom_segment(aes(x = 5.25, xend = 5.25, y = SEfitcomparisonneurologic$coefficients[1] + SEfitcomparisonneurologic$coefficients[2], yend = (SEfitcomparisonneurologic$coefficients[1] + SEfitcomparisonneurologic$coefficients[2] - SEfitcomparisonneurologic$coefficients[4])), color = "#5698a3")

## Scale for 'x' is already present. Adding another scale for 'x', which
## will replace the existing scale.

neurologichistoryfigure

# Interactive graph
ggplotly(neurologichistoryfigure)

Comparison patients with and without a history of epilepsy

## Proportion of patients with a history of epilepsy per time bin
epilepsyhistory1h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 1, ]$priorepilepsy)[2] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 1, ]$priorepilepsy)[1] + table(pSERGcircadian[pSERGcircadian$hourofSE == 1, ]$priorepilepsy)[2])
epilepsyhistory2h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 2, ]$priorepilepsy)[2] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 2, ]$priorepilepsy)[1] + table(pSERGcircadian[pSERGcircadian$hourofSE == 2, ]$priorepilepsy)[2])
epilepsyhistory3h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 3, ]$priorepilepsy)[2] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 3, ]$priorepilepsy)[1] + table(pSERGcircadian[pSERGcircadian$hourofSE == 3, ]$priorepilepsy)[2])
epilepsyhistory4h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 4, ]$priorepilepsy)[2] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 4, ]$priorepilepsy)[1] + table(pSERGcircadian[pSERGcircadian$hourofSE == 4, ]$priorepilepsy)[2])
epilepsyhistory5h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 5, ]$priorepilepsy)[2] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 5, ]$priorepilepsy)[1] + table(pSERGcircadian[pSERGcircadian$hourofSE == 5, ]$priorepilepsy)[2])
epilepsyhistory6h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 6, ]$priorepilepsy)[2] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 6, ]$priorepilepsy)[1] + table(pSERGcircadian[pSERGcircadian$hourofSE == 6, ]$priorepilepsy)[2])
epilepsyhistory7h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 7, ]$priorepilepsy)[2] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 7, ]$priorepilepsy)[1] + table(pSERGcircadian[pSERGcircadian$hourofSE == 7, ]$priorepilepsy)[2])
epilepsyhistory8h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 8, ]$priorepilepsy)[2] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 8, ]$priorepilepsy)[1] + table(pSERGcircadian[pSERGcircadian$hourofSE == 8, ]$priorepilepsy)[2])
epilepsyhistory9h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 9, ]$priorepilepsy)[2] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 9, ]$priorepilepsy)[1] + table(pSERGcircadian[pSERGcircadian$hourofSE == 9, ]$priorepilepsy)[2])
epilepsyhistory10h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 10, ]$priorepilepsy)[2] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 10, ]$priorepilepsy)[1] + table(pSERGcircadian[pSERGcircadian$hourofSE == 10, ]$priorepilepsy)[2])
epilepsyhistory11h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 11, ]$priorepilepsy)[2] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 11, ]$priorepilepsy)[1] + table(pSERGcircadian[pSERGcircadian$hourofSE == 11, ]$priorepilepsy)[2])
epilepsyhistory12h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 12, ]$priorepilepsy)[2] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 12, ]$priorepilepsy)[1] + table(pSERGcircadian[pSERGcircadian$hourofSE == 12, ]$priorepilepsy)[2])
epilepsyhistory13h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 13, ]$priorepilepsy)[2] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 13, ]$priorepilepsy)[1] + table(pSERGcircadian[pSERGcircadian$hourofSE == 13, ]$priorepilepsy)[2])
epilepsyhistory14h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 14, ]$priorepilepsy)[2] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 14, ]$priorepilepsy)[1] + table(pSERGcircadian[pSERGcircadian$hourofSE == 14, ]$priorepilepsy)[2])
epilepsyhistory15h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 15, ]$priorepilepsy)[2] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 15, ]$priorepilepsy)[1] + table(pSERGcircadian[pSERGcircadian$hourofSE == 15, ]$priorepilepsy)[2])
epilepsyhistory16h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 16, ]$priorepilepsy)[2] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 16, ]$priorepilepsy)[1] + table(pSERGcircadian[pSERGcircadian$hourofSE == 16, ]$priorepilepsy)[2])
epilepsyhistory17h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 17, ]$priorepilepsy)[2] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 17, ]$priorepilepsy)[1] + table(pSERGcircadian[pSERGcircadian$hourofSE == 17, ]$priorepilepsy)[2])
epilepsyhistory18h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 18, ]$priorepilepsy)[2] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 18, ]$priorepilepsy)[1] + table(pSERGcircadian[pSERGcircadian$hourofSE == 18, ]$priorepilepsy)[2])
epilepsyhistory19h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 19, ]$priorepilepsy)[2] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 19, ]$priorepilepsy)[1] + table(pSERGcircadian[pSERGcircadian$hourofSE == 19, ]$priorepilepsy)[2])
epilepsyhistory20h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 20, ]$priorepilepsy)[2] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 20, ]$priorepilepsy)[1] + table(pSERGcircadian[pSERGcircadian$hourofSE == 20, ]$priorepilepsy)[2])
epilepsyhistory21h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 21, ]$priorepilepsy)[2] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 21, ]$priorepilepsy)[1] + table(pSERGcircadian[pSERGcircadian$hourofSE == 21, ]$priorepilepsy)[2])
epilepsyhistory22h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 22, ]$priorepilepsy)[2] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 22, ]$priorepilepsy)[1] + table(pSERGcircadian[pSERGcircadian$hourofSE == 22, ]$priorepilepsy)[2])
epilepsyhistory23h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 23, ]$priorepilepsy)[2] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 23, ]$priorepilepsy)[1] + table(pSERGcircadian[pSERGcircadian$hourofSE == 23, ]$priorepilepsy)[2])
epilepsyhistory24h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 24, ]$priorepilepsy)[2] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 24, ]$priorepilepsy)[1] + table(pSERGcircadian[pSERGcircadian$hourofSE == 24, ]$priorepilepsy)[2])

circadiandata$epilepsy <- c(epilepsyhistory1h, epilepsyhistory2h, epilepsyhistory3h, epilepsyhistory4h, epilepsyhistory5h, epilepsyhistory6h, epilepsyhistory7h, epilepsyhistory8h, epilepsyhistory9h, epilepsyhistory10h, epilepsyhistory11h, epilepsyhistory12h, epilepsyhistory13h, epilepsyhistory14h, epilepsyhistory15h, epilepsyhistory16h, epilepsyhistory17h, epilepsyhistory18h, epilepsyhistory19h, epilepsyhistory20h, epilepsyhistory21h, epilepsyhistory22h, epilepsyhistory23h, epilepsyhistory24h)
circadiandata$epilepsy <- ifelse(is.na(circadiandata$epilepsy) == TRUE, 0, circadiandata$epilepsy)
circadiandata$numberepilepsy <- circadiandata$numberSE * circadiandata$epilepsy

## Cosinor analysis
SEfitcomparisonepilepsy <- cosinor.lm(numberSE ~ time(hourofSE) + numberepilepsy + amp.acro(numberepilepsy), data = circadiandata, period = 24)
summary(SEfitcomparisonepilepsy)

## Raw model coefficients:
##                    estimate standard.error lower.CI upper.CI p.value
## (Intercept)          6.0596         1.4519   3.2140   8.9053  0.0000
## numberepilepsy       1.0721         0.2475   0.5871   1.5571  0.0000
## rrr                 -2.1690         1.6066  -5.3179   0.9799  0.1770
## sss                  0.3465         2.6113  -4.7716   5.4646  0.8944
## numberepilepsy:rrr  -0.0495         0.2394  -0.5187   0.4197  0.8361
## numberepilepsy:sss  -0.1977         0.4082  -0.9977   0.6023  0.6281
## 
## ***********************
## 
## Transformed coefficients:
##                          estimate standard.error lower.CI upper.CI p.value
## (Intercept)                6.0596         1.4519   3.2140   8.9053  0.0000
## [numberepilepsy = 1]       1.0721         0.2475   0.5871   1.5571  0.0000
## amp                        2.1965         1.7245  -1.1834   5.5764  0.2028
## [numberepilepsy = 1]:amp   2.2235         1.4471  -0.6128   5.0599  0.1244
## acr                       -0.1584         1.1541  -2.4204   2.1036  0.8908
## [numberepilepsy = 1]:acr  -0.0670         0.9918  -2.0109   1.8770  0.9462

test_cosinor(SEfitcomparisonepilepsy, "numberepilepsy", param = "amp")

## Global test: 
## Statistic: 
## [1] 0.01
## 
## 
##  P-value: 
## [1] 0.94
## 
##  Individual tests: 
## Statistic: 
## [1] 0.08
## 
## 
##  P-value: 
## [1] 0.94
## 
##  Estimate and confidence interval[1] "0.03 (-0.68 to 0.73)"

test_cosinor(SEfitcomparisonepilepsy, "numberepilepsy", param = "acr")

## Global test: 
## Statistic: 
## [1] 0.24
## 
## 
##  P-value: 
## [1] 0.6233
## 
##  Individual tests: 
## Statistic: 
## [1] 0.49
## 
## 
##  P-value: 
## [1] 0.6233
## 
##  Estimate and confidence interval[1] "0.09 (-0.27 to 0.46)"

##Graph
epilepsyhistoryfigure <- ggplot.cosinor.lm(SEfitcomparisonepilepsy, x_str = "numberepilepsy") +
  geom_bar(data = pSERGcircadian[pSERGcircadian$priorepilepsy == 1, ], aes(x = hourofSE, y = events), stat = 'identity', color = "#5698a3", fill = "#5698a3", alpha = 0.2) + 
  geom_bar(data = pSERGcircadian[pSERGcircadian$priorepilepsy == 0, ], aes(x = hourofSE, y = events), stat = 'identity', color = "#a30234", fill = "#a30234", alpha = 0.2) +
    scale_x_discrete(limits = c(1,24)) +
  theme(panel.background=element_rect(fill = "white"), 
        axis.text=element_text(size = 18, color = "black", face = "bold"),
        axis.title=element_text(size = 24, color = "black", face = "bold"),
        axis.ticks.length=unit(0.4, "cm"), axis.ticks = element_line(size = 1.2),
        legend.position = "none") +
  scale_x_continuous(breaks = c(4, 8, 12, 16, 20, 24)) +
  scale_y_continuous(breaks = c(5, 10, 15)) +
  labs(x = "Time of day (hour)", y = "Number of rSE episodes") +
  geom_segment(aes(x = 0, xend = 24, y = SEfitcomparisonepilepsy$coefficients[1], yend = SEfitcomparisonepilepsy$coefficients[1]), color = "#a30234") +
  geom_segment(aes(x = 11.35, xend = 11.35, y = SEfitcomparisonepilepsy$coefficients[1], yend = (SEfitcomparisonepilepsy$coefficients[1] + SEfitcomparisonepilepsy$coefficients[3])), color = "#a30234") +
  geom_segment(aes(x = 23.35, xend = 23.35, y = SEfitcomparisonepilepsy$coefficients[1], yend = (SEfitcomparisonepilepsy$coefficients[1] - SEfitcomparisonepilepsy$coefficients[3])), color = "#a30234") +
  geom_segment(aes(x = 0, xend = 24, y = SEfitcomparisonepilepsy$coefficients[1] + SEfitcomparisonepilepsy$coefficients[2], yend = SEfitcomparisonepilepsy$coefficients[1] + SEfitcomparisonepilepsy$coefficients[2]), color = "#5698a3") +
  geom_segment(aes(x = 11.7, xend = 11.7, y = SEfitcomparisonepilepsy$coefficients[1] + SEfitcomparisonepilepsy$coefficients[2], yend = (SEfitcomparisonepilepsy$coefficients[1] + SEfitcomparisonepilepsy$coefficients[2] + SEfitcomparisonepilepsy$coefficients[4])), color = "#5698a3") + geom_segment(aes(x = 23.7, xend = 23.7, y = SEfitcomparisonepilepsy$coefficients[1] + SEfitcomparisonepilepsy$coefficients[2], yend = (SEfitcomparisonepilepsy$coefficients[1] + SEfitcomparisonepilepsy$coefficients[2] - SEfitcomparisonepilepsy$coefficients[4])), color = "#5698a3")

## Scale for 'x' is already present. Adding another scale for 'x', which
## will replace the existing scale.

epilepsyhistoryfigure

# Interactive graph
ggplotly(epilepsyhistoryfigure)

Comparison patients with SE onset in and out of the hospital

## Proportion of patients with onset of SE out of the hospital
outofhospitalonset1h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 1, ]$HOSPITALONSET)[2] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 1, ]$HOSPITALONSET)[3] + table(pSERGcircadian[pSERGcircadian$hourofSE == 1, ]$HOSPITALONSET)[2])
outofhospitalonset2h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 2, ]$HOSPITALONSET)[2] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 2, ]$HOSPITALONSET)[3] + table(pSERGcircadian[pSERGcircadian$hourofSE == 2, ]$HOSPITALONSET)[2])
outofhospitalonset3h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 3, ]$HOSPITALONSET)[2] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 3, ]$HOSPITALONSET)[3] + table(pSERGcircadian[pSERGcircadian$hourofSE == 3, ]$HOSPITALONSET)[2])
outofhospitalonset4h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 4, ]$HOSPITALONSET)[2] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 4, ]$HOSPITALONSET)[3] + table(pSERGcircadian[pSERGcircadian$hourofSE == 4, ]$HOSPITALONSET)[2])
outofhospitalonset5h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 5, ]$HOSPITALONSET)[2] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 5, ]$HOSPITALONSET)[3] + table(pSERGcircadian[pSERGcircadian$hourofSE == 5, ]$HOSPITALONSET)[2])
outofhospitalonset6h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 6, ]$HOSPITALONSET)[2] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 6, ]$HOSPITALONSET)[3] + table(pSERGcircadian[pSERGcircadian$hourofSE == 6, ]$HOSPITALONSET)[2])
outofhospitalonset7h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 7, ]$HOSPITALONSET)[2] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 7, ]$HOSPITALONSET)[3] + table(pSERGcircadian[pSERGcircadian$hourofSE == 7, ]$HOSPITALONSET)[2])
outofhospitalonset8h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 8, ]$HOSPITALONSET)[2] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 8, ]$HOSPITALONSET)[3] + table(pSERGcircadian[pSERGcircadian$hourofSE == 8, ]$HOSPITALONSET)[2])
outofhospitalonset9h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 9, ]$HOSPITALONSET)[2] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 9, ]$HOSPITALONSET)[3] + table(pSERGcircadian[pSERGcircadian$hourofSE == 9, ]$HOSPITALONSET)[2])
outofhospitalonset10h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 10, ]$HOSPITALONSET)[2] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 10, ]$HOSPITALONSET)[3] + table(pSERGcircadian[pSERGcircadian$hourofSE == 10, ]$HOSPITALONSET)[2])
outofhospitalonset11h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 11, ]$HOSPITALONSET)[2] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 11, ]$HOSPITALONSET)[3] + table(pSERGcircadian[pSERGcircadian$hourofSE == 11, ]$HOSPITALONSET)[2])
outofhospitalonset12h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 12, ]$HOSPITALONSET)[2] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 12, ]$HOSPITALONSET)[3] + table(pSERGcircadian[pSERGcircadian$hourofSE == 12, ]$HOSPITALONSET)[2])
outofhospitalonset13h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 13, ]$HOSPITALONSET)[2] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 13, ]$HOSPITALONSET)[3] + table(pSERGcircadian[pSERGcircadian$hourofSE == 13, ]$HOSPITALONSET)[2])
outofhospitalonset14h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 14, ]$HOSPITALONSET)[2] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 14, ]$HOSPITALONSET)[3] + table(pSERGcircadian[pSERGcircadian$hourofSE == 14, ]$HOSPITALONSET)[2])
outofhospitalonset15h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 15, ]$HOSPITALONSET)[2] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 15, ]$HOSPITALONSET)[3] + table(pSERGcircadian[pSERGcircadian$hourofSE == 15, ]$HOSPITALONSET)[2])
outofhospitalonset16h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 16, ]$HOSPITALONSET)[2] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 16, ]$HOSPITALONSET)[3] + table(pSERGcircadian[pSERGcircadian$hourofSE == 16, ]$HOSPITALONSET)[2])
outofhospitalonset17h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 17, ]$HOSPITALONSET)[2] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 17, ]$HOSPITALONSET)[3] + table(pSERGcircadian[pSERGcircadian$hourofSE == 17, ]$HOSPITALONSET)[2])
outofhospitalonset18h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 18, ]$HOSPITALONSET)[2] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 18, ]$HOSPITALONSET)[3] + table(pSERGcircadian[pSERGcircadian$hourofSE == 18, ]$HOSPITALONSET)[2])
outofhospitalonset19h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 19, ]$HOSPITALONSET)[2] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 19, ]$HOSPITALONSET)[3] + table(pSERGcircadian[pSERGcircadian$hourofSE == 19, ]$HOSPITALONSET)[2])
outofhospitalonset20h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 20, ]$HOSPITALONSET)[2] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 20, ]$HOSPITALONSET)[3] + table(pSERGcircadian[pSERGcircadian$hourofSE == 20, ]$HOSPITALONSET)[2])
outofhospitalonset21h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 21, ]$HOSPITALONSET)[2] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 21, ]$HOSPITALONSET)[3] + table(pSERGcircadian[pSERGcircadian$hourofSE == 21, ]$HOSPITALONSET)[2])
outofhospitalonset22h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 22, ]$HOSPITALONSET)[2] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 22, ]$HOSPITALONSET)[3] + table(pSERGcircadian[pSERGcircadian$hourofSE == 22, ]$HOSPITALONSET)[2])
outofhospitalonset23h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 23, ]$HOSPITALONSET)[2] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 23, ]$HOSPITALONSET)[3] + table(pSERGcircadian[pSERGcircadian$hourofSE == 23, ]$HOSPITALONSET)[2])
outofhospitalonset24h <- table(pSERGcircadian[pSERGcircadian$hourofSE == 24, ]$HOSPITALONSET)[2] / (table(pSERGcircadian[pSERGcircadian$hourofSE == 24, ]$HOSPITALONSET)[3] + table(pSERGcircadian[pSERGcircadian$hourofSE == 24, ]$HOSPITALONSET)[2])


circadiandata$outofhospitalonset <- c(outofhospitalonset1h, outofhospitalonset2h, outofhospitalonset3h,
                                      outofhospitalonset4h, outofhospitalonset5h, outofhospitalonset6h,
                                      outofhospitalonset7h, outofhospitalonset8h, outofhospitalonset9h,
                                      outofhospitalonset10h, outofhospitalonset11h, outofhospitalonset12h,
                                      outofhospitalonset13h, outofhospitalonset14h, outofhospitalonset15h,
                                      outofhospitalonset16h, outofhospitalonset17h, outofhospitalonset18h,
                                      outofhospitalonset19h, outofhospitalonset20h, outofhospitalonset21h,
                                      outofhospitalonset22h, outofhospitalonset23h, outofhospitalonset24h)
circadiandata$outofhospitalonset <- ifelse(is.na(circadiandata$outofhospitalonset) == TRUE, 0, circadiandata$outofhospitalonset)
circadiandata$numberoutofhospitalonset <- circadiandata$numberSE * circadiandata$outofhospitalonset

## Cosinor analysis
SEfitcomparisonhospitalonset <- cosinor.lm(numberSE ~ time(hourofSE) + numberoutofhospitalonset + amp.acro(numberoutofhospitalonset), data = circadiandata, period = 24)
summary(SEfitcomparisonhospitalonset)

## Raw model coefficients:
##                              estimate standard.error lower.CI upper.CI
## (Intercept)                    4.4437         1.9984   0.5269   8.3606
## numberoutofhospitalonset       0.9528         0.2308   0.5005   1.4051
## rrr                           -0.8620         2.1684  -5.1120   3.3880
## sss                           -0.1296         3.5156  -7.0202   6.7609
## numberoutofhospitalonset:rrr  -0.0420         0.2561  -0.5438   0.4599
## numberoutofhospitalonset:sss   0.0637         0.3804  -0.6820   0.8093
##                              p.value
## (Intercept)                   0.0262
## numberoutofhospitalonset      0.0000
## rrr                           0.6910
## sss                           0.9706
## numberoutofhospitalonset:rrr  0.8698
## numberoutofhospitalonset:sss  0.8671
## 
## ***********************
## 
## Transformed coefficients:
##                                    estimate standard.error lower.CI
## (Intercept)                          4.4437         1.9984   0.5269
## [numberoutofhospitalonset = 1]       0.9528         0.2308   0.5005
## amp                                  0.8717         2.1778  -3.3967
## [numberoutofhospitalonset = 1]:amp   0.9064         1.9308  -2.8778
## acr                                  0.1493         4.0264  -7.7423
## [numberoutofhospitalonset = 1]:acr   0.0728         3.4760  -6.7400
##                                    upper.CI p.value
## (Intercept)                          8.3606  0.0262
## [numberoutofhospitalonset = 1]       1.4051  0.0000
## amp                                  5.1401  0.6890
## [numberoutofhospitalonset = 1]:amp   4.6906  0.6387
## acr                                  8.0408  0.9704
## [numberoutofhospitalonset = 1]:acr   6.8856  0.9833

test_cosinor(SEfitcomparisonhospitalonset, "numberoutofhospitalonset", param = "amp")

## Global test: 
## Statistic: 
## [1] 0.01
## 
## 
##  P-value: 
## [1] 0.9248
## 
##  Individual tests: 
## Statistic: 
## [1] 0.09
## 
## 
##  P-value: 
## [1] 0.9248
## 
##  Estimate and confidence interval[1] "0.03 (-0.69 to 0.76)"

test_cosinor(SEfitcomparisonhospitalonset, "numberoutofhospitalonset", param = "acr")

## Global test: 
## Statistic: 
## [1] 0.02
## 
## 
##  P-value: 
## [1] 0.896
## 
##  Individual tests: 
## Statistic: 
## [1] -0.13
## 
## 
##  P-value: 
## [1] 0.896
## 
##  Estimate and confidence interval[1] "-0.08 (-1.22 to 1.07)"

##Graph
hospitalonsetfigure <- ggplot.cosinor.lm(SEfitcomparisonhospitalonset, x_str = "numberoutofhospitalonset") +
  geom_bar(data = pSERGcircadian[pSERGcircadian$HOSPITALONSET == "no", ], aes(x = hourofSE, y = events), stat = 'identity', color = "#5698a3", fill = "#5698a3", alpha = 0.2) + 
  geom_bar(data = pSERGcircadian[pSERGcircadian$HOSPITALONSET == "yes", ], aes(x = hourofSE, y = events), stat = 'identity', color = "#a30234", fill = "#a30234", alpha = 0.2) +
    scale_x_discrete(limits = c(1,24)) +
  theme(panel.background=element_rect(fill = "white"), 
        axis.text=element_text(size = 18, color = "black", face = "bold"),
        axis.title=element_text(size = 24, color = "black", face = "bold"),
        axis.ticks.length=unit(0.4, "cm"), axis.ticks = element_line(size = 1.2),
        legend.position = "none") +
  scale_x_continuous(breaks = c(4, 8, 12, 16, 20, 24)) +
  scale_y_continuous(breaks = c(5, 10, 15)) +
  labs(x = "Time of day (hour)", y = "Number of rSE episodes") +
  geom_segment(aes(x = 0, xend = 24, y = SEfitcomparisonhospitalonset$coefficients[1], yend = SEfitcomparisonhospitalonset$coefficients[1]), color = "#a30234") +
  geom_segment(aes(x = 12.5, xend = 12.5, y = SEfitcomparisonhospitalonset$coefficients[1], yend = (SEfitcomparisonhospitalonset$coefficients[1] + SEfitcomparisonhospitalonset$coefficients[3])), color = "#a30234") +
  geom_segment(aes(x = 0.5, xend = 0.5, y = SEfitcomparisonhospitalonset$coefficients[1], yend = (SEfitcomparisonhospitalonset$coefficients[1] - SEfitcomparisonhospitalonset$coefficients[3])), color = "#a30234") +
  geom_segment(aes(x = 0, xend = 24, y = SEfitcomparisonhospitalonset$coefficients[1] + SEfitcomparisonhospitalonset$coefficients[2], yend = SEfitcomparisonhospitalonset$coefficients[1] + SEfitcomparisonhospitalonset$coefficients[2]), color = "#5698a3") +
  geom_segment(aes(x = 12.2, xend = 12.2, y = SEfitcomparisonhospitalonset$coefficients[1] + SEfitcomparisonhospitalonset$coefficients[2], yend = (SEfitcomparisonhospitalonset$coefficients[1] + SEfitcomparisonhospitalonset$coefficients[2] + SEfitcomparisonhospitalonset$coefficients[4])), color = "#5698a3") + geom_segment(aes(x = 0.2, xend = 0.2, y = SEfitcomparisonhospitalonset$coefficients[1] + SEfitcomparisonhospitalonset$coefficients[2], yend = (SEfitcomparisonhospitalonset$coefficients[1] + SEfitcomparisonhospitalonset$coefficients[2] - SEfitcomparisonhospitalonset$coefficients[4])), color = "#5698a3")

## Scale for 'x' is already present. Adding another scale for 'x', which
## will replace the existing scale.

hospitalonsetfigure

# Interactive graph
ggplotly(hospitalonsetfigure)

Time to first BZD

##Transform hours of SE to factor
pSERGcircadian$hourofSE<- as.factor(pSERGcircadian$hourofSE)

# Transform time to first BZD time to numeric
pSERGcircadian$BZDTIME.0 <- as.numeric(as.character(pSERGcircadian$BZDTIME.0))

## Warning: NAs introduced by coercion

#Summary data
summary(pSERGcircadian$BZDTIME.0)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##    0.00    5.00   17.50   74.58   49.75 1440.00      22

######################Transform hours of SE to numeric############
pSERGcircadian$hourofSE <- as.numeric(pSERGcircadian$hourofSE)

#Cosinor analysis
BZDfit <- cosinor.lm(BZDTIME.0 ~ time(hourofSE), data = pSERGcircadian[!is.na(pSERGcircadian$hourofSE), ], period = 24)
summary(BZDfit)

## Raw model coefficients:
##             estimate standard.error lower.CI upper.CI p.value
## (Intercept)  75.4565        11.7140  52.4975  98.4156  0.0000
## rrr           9.7029        16.7027 -23.0338  42.4396  0.5613
## sss           1.7495        16.3375 -30.2714  33.7705  0.9147
## 
## ***********************
## 
## Transformed coefficients:
##             estimate standard.error lower.CI upper.CI p.value
## (Intercept)  75.4565        11.7140  52.4975  98.4156  0.0000
## amp           9.8594        16.8564 -23.1786  42.8973  0.5586
## acr           0.1784         1.6410  -3.0378   3.3946  0.9134

##Draw the number of SE in a plot
################Transform hours of SE to factor######
pSERGcircadian$hourofSE <- as.factor(pSERGcircadian$hourofSE)
timeBZD <- ggplot.cosinor.lm(BZDfit) +
  geom_boxplot(data = pSERGcircadian[!is.na(pSERGcircadian$hourofSE), ], aes(x = hourofSE, y = BZDTIME.0), color = "#0076c0", fill = "#67771a", alpha = 0.5) + 
  scale_x_discrete(limits = c(1,24)) + 
  theme(panel.background=element_rect(fill = "white"), 
        axis.text=element_text(size = 16, color = "black", face = "bold"),
        axis.title=element_text(size = 18, color = "black", face = "bold"),
        axis.ticks.length = unit(0.4, "cm"), axis.ticks = element_line(size = 1.2)) +
  scale_x_discrete(breaks = c("4", "8", "12", "16", "20", "24")) + 
  scale_y_continuous(breaks = c(0, 50, 100, 150, 200, 250, 300, 350)) + 
  coord_cartesian(ylim = c(0, 200)) +
  labs(x = "Time of day (hour)", y= "Time to the first BZD (min)") +
  geom_segment(aes(x = 0,xend = 24, y= BZDfit$coefficients[1], yend = BZDfit$coefficients[1])) +
  geom_segment(aes(x = 0.75, xend = 0.75, y = BZDfit$coefficients[1], yend = BZDfit$coefficients[1] + BZDfit$coefficients[2])) +
  geom_segment(aes(x = 12.75, xend = 12.75, y = BZDfit$coefficients[1], yend = BZDfit$coefficients[1] - BZDfit$coefficients[2]))

## Scale for 'x' is already present. Adding another scale for 'x', which
## will replace the existing scale.

timeBZD

## Warning: Removed 22 rows containing non-finite values (stat_boxplot).

# Interactive graph
ggplotly(timeBZD)

## Warning: Removed 22 rows containing non-finite values (stat_boxplot).

Time to first non-BZD AED

## Transform hours of SE to factor
pSERGcircadian$hourofSE <- as.factor(pSERGcircadian$hourofSE)


# Transform time to first AED time to numeric
pSERGcircadian$AEDTIME.0 <- as.numeric(as.character(pSERGcircadian$AEDTIME.0))

## Warning: NAs introduced by coercion

# Summary data
summary(pSERGcircadian$AEDTIME.0)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##     0.0    35.0    66.0   174.7   155.0  4320.0      23

## Transform hours of SE to numeric 
pSERGcircadian$hourofSE <- as.numeric(pSERGcircadian$hourofSE)

# Cosinor analysis
ASMfit <- cosinor.lm(AEDTIME.0~time(hourofSE), data=pSERGcircadian, period=24)
summary(ASMfit)

## Raw model coefficients:
##             estimate standard.error lower.CI upper.CI p.value
## (Intercept) 178.8390        23.2214 133.3258 224.3521  0.0000
## rrr          45.2087        33.5534 -20.5547 110.9721  0.1779
## sss          56.2470        31.9973  -6.4666 118.9606  0.0788
## 
## ***********************
## 
## Transformed coefficients:
##             estimate standard.error lower.CI upper.CI p.value
## (Intercept) 178.8390        23.2214 133.3258 224.3521  0.0000
## amp          72.1634        33.7732   5.9692 138.3576  0.0326
## acr           0.8938         0.4402   0.0310   1.7565  0.0423

################Transform hours of SE to factor######
pSERGcircadian$hourofSE <- as.factor(pSERGcircadian$hourofSE)

#Draw the data and the fitting curve
timeASM <- ggplot.cosinor.lm(ASMfit) +
  geom_boxplot(data = pSERGcircadian[!is.na(pSERGcircadian$hourofSE), ], aes(hourofSE, AEDTIME.0), color = "#0076c0", fill = "#67771a", alpha = 0.5) +
  scale_x_discrete(limits = c(1,24)) + 
  scale_y_discrete(limits = c(0,800)) + 
  theme(panel.background = element_rect(fill = "white"),
        axis.text = element_text(size = 16, color = "black", face = "bold"),
        axis.title = element_text(size = 18, color = "black", face = "bold"),
        axis.ticks.length = unit(0.4, "cm"), axis.ticks = element_line(size = 1.2)) +
  scale_x_discrete(breaks = c("4", "8", "12", "16", "20", "24")) +
  scale_y_continuous(breaks = c(0, 100, 200, 300, 400, 500, 600)) + 
  coord_cartesian(ylim = c(0, 600)) +
  labs(x = "Time of day (hour)", y = "Time to the first non-BZD AED (min)") +
  geom_segment(aes(x = 0, xend = 24, y = ASMfit$coefficients[1], yend = ASMfit$coefficients[1])) +
  geom_segment(aes(x = 3.3, xend = 3.3, y = ASMfit$coefficients[1], yend = ASMfit$coefficients[1] + ASMfit$coefficients[2])) +
  geom_segment(aes(x = 15.3, xend = 15.3, y = ASMfit$coefficients[1], yend = ASMfit$coefficients[1] - ASMfit$coefficients[2]))

## Scale for 'x' is already present. Adding another scale for 'x', which
## will replace the existing scale.

## Scale for 'y' is already present. Adding another scale for 'y', which
## will replace the existing scale.

timeASM

## Warning: Removed 23 rows containing non-finite values (stat_boxplot).

# Interactive graph
ggplotly(timeASM)

## Warning: Removed 23 rows containing non-finite values (stat_boxplot).

Alternative hour definition

#################Create variable sensitivity analysis hour of SE######################
## Define the newly defined hour
pSERGcircadian$minutes <- pSERGcircadian$TIMESEIZURE_MINUTES
pSERGcircadian$minutesashours <- pSERGcircadian$minutes / 60
pSERGcircadian$minutesashours <- ifelse(is.na(pSERGcircadian$minutesashours) == TRUE | pSERGcircadian$minutesashours > 1, 0, pSERGcircadian$minutesashours)

## Create variable new hour
pSERGcircadian$newhour <- as.numeric(pSERGcircadian$hourofSE) + pSERGcircadian$minutesashours

## Rename hours
pSERGcircadian$newhour[pSERGcircadian$newhour >= 24.5 & pSERGcircadian$newhour < 1.5] <- 1
pSERGcircadian$newhour[pSERGcircadian$newhour >= 1.5 & pSERGcircadian$newhour < 2.5] <- 2
pSERGcircadian$newhour[pSERGcircadian$newhour >= 2.5 & pSERGcircadian$newhour < 3.5] <- 3
pSERGcircadian$newhour[pSERGcircadian$newhour >= 3.5 & pSERGcircadian$newhour < 4.5] <- 4
pSERGcircadian$newhour[pSERGcircadian$newhour >= 4.5 & pSERGcircadian$newhour < 5.5] <- 5
pSERGcircadian$newhour[pSERGcircadian$newhour >= 5.5 & pSERGcircadian$newhour < 6.5] <- 6
pSERGcircadian$newhour[pSERGcircadian$newhour >= 6.5 & pSERGcircadian$newhour < 7.5] <- 7
pSERGcircadian$newhour[pSERGcircadian$newhour >= 7.5 & pSERGcircadian$newhour < 8.5] <- 8
pSERGcircadian$newhour[pSERGcircadian$newhour >= 8.5 & pSERGcircadian$newhour < 9.5] <- 9
pSERGcircadian$newhour[pSERGcircadian$newhour >= 9.5 & pSERGcircadian$newhour < 10.5] <- 10
pSERGcircadian$newhour[pSERGcircadian$newhour >= 10.5 & pSERGcircadian$newhour < 11.5] <- 11
pSERGcircadian$newhour[pSERGcircadian$newhour >= 11.5 & pSERGcircadian$newhour < 12.5] <- 12
pSERGcircadian$newhour[pSERGcircadian$newhour >= 12.5 & pSERGcircadian$newhour < 13.5] <- 13
pSERGcircadian$newhour[pSERGcircadian$newhour >= 13.5 & pSERGcircadian$newhour < 14.5] <- 14
pSERGcircadian$newhour[pSERGcircadian$newhour >= 14.5 & pSERGcircadian$newhour < 15.5] <- 15
pSERGcircadian$newhour[pSERGcircadian$newhour >= 15.5 & pSERGcircadian$newhour < 16.5] <- 16
pSERGcircadian$newhour[pSERGcircadian$newhour >= 16.5 & pSERGcircadian$newhour < 17.5] <- 17
pSERGcircadian$newhour[pSERGcircadian$newhour >= 17.5 & pSERGcircadian$newhour < 18.5] <- 18
pSERGcircadian$newhour[pSERGcircadian$newhour >= 18.5 & pSERGcircadian$newhour < 19.5] <- 19
pSERGcircadian$newhour[pSERGcircadian$newhour >= 19.5 & pSERGcircadian$newhour < 20.5] <- 19
pSERGcircadian$newhour[pSERGcircadian$newhour >= 20.5 & pSERGcircadian$newhour < 21.5] <- 20
pSERGcircadian$newhour[pSERGcircadian$newhour >= 21.5 & pSERGcircadian$newhour < 22.5] <- 21
pSERGcircadian$newhour[pSERGcircadian$newhour >= 22.5 & pSERGcircadian$newhour < 23.5] <- 22
pSERGcircadian$newhour[pSERGcircadian$newhour >= 23.5 & pSERGcircadian$newhour < 24.5] <- 24


# Test (Kolmogorov-Smirnov) if the distribution of SE during the day is consistent with a uniform distribution
ks.test(pSERGcircadian$newhour, "punif", 0, 24)

## Warning in ks.test(pSERGcircadian$newhour, "punif", 0, 24): ties should not
## be present for the Kolmogorov-Smirnov test

## 
##  One-sample Kolmogorov-Smirnov test
## 
## data:  pSERGcircadian$newhour
## D = 0.10333, p-value = 0.003301
## alternative hypothesis: two-sided

#####################################################################################
# Create another database with the circadian data
newoccurrences <- count(pSERGcircadian, "newhour")
# Make it a dataframe
newcircadiandata <- data.frame(newoccurrences)
newcircadiandata$numberSE <- newcircadiandata$freq
newcircadiandata <- newcircadiandata[ , c("newhour", "numberSE")]
#####################################################################################



# Cosinor analysis
newSEfit <- cosinor.lm(numberSE ~ time(newhour), data = newcircadiandata, period = 24)
summary(newSEfit)

## Raw model coefficients:
##             estimate standard.error lower.CI upper.CI p.value
## (Intercept)  11.1711         0.9036   9.4001  12.9422  0.0000
## rrr          -5.7450         1.1758  -8.0495  -3.4405  0.0000
## sss          -0.7859         1.3776  -3.4859   1.9141  0.5683
## 
## ***********************
## 
## Transformed coefficients:
##             estimate standard.error lower.CI upper.CI p.value
## (Intercept)  11.1711         0.9036   9.4001  12.9422  0.0000
## amp           5.7985         1.1635   3.5182   8.0789  0.0000
## acr           0.1360         0.2394  -0.3332   0.6051  0.5701

################Every hour new definition ###############
## Draw the number of SE episodes in a plot with stacked units
newdistributionhour <- ggplot(pSERGcircadian, aes(x = round(newhour), y = events)) + geom_bar(stat = 'identity', color = "#0076c0", fill = "#67771a") + 
  scale_x_discrete(limits = c(1,24)) +
  coord_cartesian(xlim = c(1, 23.4)) +
  theme(panel.background = element_rect(fill = "white"), 
        axis.text = element_text(size = 18, color = "black", face = "bold"),
        axis.title = element_text(size = 24, color = "black", face = "bold"),
        axis.ticks.length = unit(0.4, "cm"), axis.ticks = element_line(size = 1.2)) +
  scale_x_continuous(breaks = c(4, 8, 12, 16, 20, 24)) +
  scale_y_continuous(breaks = c(5, 10, 15, 20, 25)) +
  labs(x= "Time of day (hour)", y= "Number of rSE episodes")

## Scale for 'x' is already present. Adding another scale for 'x', which
## will replace the existing scale.

newdistributionhour

# Interactive graph
ggplotly(newdistributionhour)